Dynamics of a Massive Black Hole at the Center of a Dense 

Stellar System 

Pinaki Chatterjee^, Lars Hernquist^ & Abraham Loeb^ 
Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 0213^ 

§ ■ ABSTRACT 

(N 

<D • We develop a simple physical model to describe the dynamics of a massive 

'■^ ■ point-like object, such as a black hole, near the center of a dense stellar 

system. It is shown that the total force on this body can be separated into two 
independent parts, one of which is the slowly varying influence of the aggregate 

^ ■ stellar system, and the other being the rapidly fluctuating stochastic force due to 

CN ' system distributed according to a Plummer model, it is shown that the motion 



o 



of the black hole is then similar to that of a Brownian particle in a harmonic 
potential, and we analyze its dynamics using an approach akin to Langevin's 
solution of the Brownian motion problem. The equations are solved to obtain 
Q_|. the average values, time-autocorrelation functions, and probability distributions 

Q . of the black hole's position and velocity. By comparing these results to N-body 

^ . simulations, we demonstrate that this model provides a very good statistical 

d • description of the actual black hole dynamics. As an application of our model. 



we use our results to derive a lower limit on the mass of the black hole Sgr A 



* 



^^ . in the Galactic center. 
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1. Introduction 



Black holes are thought to be ubiquitous in dense stellar systems. Matter accreting 
onto supermassive black holes near the centers of galaxies is believed to be responsible for 
the energetic emission produced by active galactic nuclei (Zel'dovich 1964; Salpeter 1964; 
Lynden-Bell 1969; Rees 1984). Furthermore, it has been conjectured that all galaxies harbor 
such black holes at their centers (but see Gebhardt et al. 2001 for recent observations 
that some do not). Although definitive proof of this hypothesis is still lacking, there exists 
evidence in some galaxies, such as NGC 4258 (Greenhill et al. 1995; Kormendy & Richstone 
1995) and our own Galaxy (see Melia & Falcke 2001 for a review), for the presence of an 
unresolved central dark mass of such high density that it is unlikely to be anything other 
than a black hole (Maoz 1998). In the case of the Galactic center, which is thought to 
coincide with the unusual radio source Sgr A*, future observations will measure the orbits of 
individual stars within 0.1" of Sgr A* (see, e.g., Ghez et al. 2000). In addition, forthcoming 
radio observations will significantly improve the current limits on the proper motion of Sgr 
A* itself (M. Reid 2001, private communication). It is important, therefore, to understand 
the general properties of the dynamics of massive bodies in dense stellar systems so that the 
observations can be unambiguously interpreted and predictions can be made to stringently 
test underlying theories. 

To pursue this goal, we present a simple model for the dynamics of a single massive 
black hole at the center of a dense stellar system. Our approach is motivated by the 
recognition (Chandrasekhar 1943a) that the force acting on an object in a stellar system 
broadly consists of two independent contributions: one part, which originates from the 
"smoothed-out" average distribution of matter in the stellar system, will vary slowly with 
position and time; the second part, which arises from discrete encounters with individual 
stars, will fluctuate much more rapidly. 

The smooth force itself is expected to be made up of two pieces: the first is the force 
arising from the potential of the aggregate distribution of stars at the position of the object; 
and the second is the dissipative force known as dynamical friction, which causes the object 
to decelerate as it moves through the stellar background (Chandrasekhar 1943b). 

The problem of the dynamics of a black hole in a stellar system is then similar 
in spirit to the Langevin model of Brownian motion (see, e.g., Chandrasekhar 1943a), 
which describes the irregular motions suffered by dust grains immersed in a gas. In the 
Langevin analysis, a Brownian particle experiences a decelerating force due to friction 
which is proportional to its velocity, and it experiences an essentially random, rapidly 
fluctuating force owing to the large rate of collisions it suffers with the gas molecules in its 
neighborhood. 
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We extend this method of analysis to the black hole problem. We take the stellar 
system to be distributed according to a Plummer potential (see Binney & Tremaine 1987, 
hereafter BT) because the dynamical equations are then relatively tractable, and because 
this density profile provides a reasonably good fit to actual stellar systems. In §2, we set up 
the model equations and also provide a justification for breaking up the force on the black 
hole into two independent parts, one smooth and slowly varying, and the other rapidly 
fiuctuating. The equations of motion for the black hole are shown to be similar to those of a 
Brownian particle in a harmonic potential well. In §3, we solve the equations of motion for 
the average position and velocity of the black hole, and the time-autocorrelation function of 
its position and velocity, obtaining both the transient and steady-state components of these 
functions. In §4, we derive the probability distributions of the black hole's position and 
velocity by solving the Fokker-Planck equation of the model. It is shown that in the steady 
state, these two variables are distributed independently with a Gaussian distribution. The 
conclusions of §3 and §4 are tested in §5 by comparing them with the results of N-body 
simulations of various systems. In §6, we combine the results of the model and observational 
limits on the proper motion of Sgr A* with physical arguments relating to the maximum 
lifetime of the cluster of stars surrounding the black hole (following the approach in Maoz 
1998) to derive lower limits on the mass of Sgr A*. Finally, §7 summarizes the paper. 



2. The Model 

Consider a black hole of mass m in a cluster of stars which we take to be described by 
a Plummer model of total mass M and length parameter a. Thus, the density and potential 
profiles are given, respectively, by 

. X _ 3Ma^ 1 

where G is the gravitational constant and r is the radial position vector from the center of 
the stellar system, which is taken as the origin. The total mass in stars inside radius r is 
then 

^"ir) = i;j^^^. (3) 

Given the potential and density profiles, one can calculate the phase space distribution 
function /, which in general depends both on position r, and velocity v, and which is 
defined such that /(r, v) d^rd^v is the mass in stars in the phase space volume d^rd^v. We 
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make the assumption that for the spherically symmetric Plummer model, / is a function of 
the relative energy per unit mass E only (and independent of specific angular momentum), 
where E- = —\v'^ — $(r) = ^(r) — ^f^, \l'(r) = — $(r) being the relative potential. The 
distribution function can then be calculated by the following equation (see BT): 



m ' '""'' "'^ 



For the Plummer case, we get 

With these preliminaries, we are in a position to calculate the forces on the black hole 
in this model. There are three such forces: the restoring force of the stellar potential, 
dynamical friction, and a random force due to discrete encounters with stars. 

The restoring force on the black hole of the stellar potential is given by F = — mV$(r), 
where r is the position vector of the black hole. Now 

GM GM f r2 3r^ 



(r2 + a2)i/2 a V 2a^ 



ia 



Since the black hole is much more massive than the stars, its typical excursion from the 
center r is small compared with a, and we are entitled to neglect terms in the above 
equation of higher order than r^. The dominant restoring force on the black hole thus takes 
the form of Hooke's law: 

F = -kr, (5) 

where the "spring constant" k is given by 

k=^. (6) 

As the black hole moves through the sea of stars, it experiences a force of deceleration 
known as dynamical friction. We use for this the Chandrasekhar dynamical friction formula 
(Chandrasekhar 1943b, BT): 

F = -/?v. 



where 



/Q ifi 21 A r^2 / , ^lo fir,u)u'^du 
(j = loTT InA Gm[m + m^j — 



y3 



In the above, v is the velocity of the black hole, m^ is the mass of each star (in the following, 
we take all stars to have equal masses, for simplicity), and InA is the Coulomb logarithm, 
which will be calculated below. 
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Note that the above formula was originally derived for the case of a mass moving 
through a homogeneous stellar system, for which the distribution function would be 
independent of r. However, it is a good approximation to replace this in the case of 
non-homogeneous systems with the distribution function in the vicinity of the black hole 
(see BT), especially since the distribution function for the Plummer model varies slowly 
with r in the region in which the black hole hole is confined {r <^ a). Since the black hole 
moves very slowly compared with the stars, we may replace /(r, u) in the integral by /(r, 0) 
to obtain (see BT): 

(3 = InA G'^mim + m^)f{r, 0). 

But 

96 Ma^ , ,, ,,7/9 

^'''•" = °' = 778;5(GmF<-*<''" ■ 
and for r <^ a^ ^(r^ ':^ —GM/a. Thus, we finally get 

/3=-^lnA(^^ m(m + m.), 



or, since m 3> rrii,, 



^ 128v^, A / G- V/2 2 

'' = ^^'"^(m?) "■ '"> 



The factor A in the Coulomb logarithm is given by 



r2 

^-"h~-"~Gm ' 



h h V^ 



where bmax and bmin are, respectively, the maximum and minimum impact parameters 
between the black hole and the stars that need be considered; bmin is usually set to be (see 
BT, Maoz 1993) bmin ~ Gitl/Vq, where Vq is the typical relative speed between the black 
hole and the stars with which it interacts. Since the velocity of the black hole is much 
smaller than that of the stars, we set Vq to be the mean squared velocity of the stars: 

Jv^f{S)47rv^dv ^ .^ GM 



The maximum impact parameter bmax is not well-defined; however, an error in the choice of 
bmax results in a much smaller error in the coefficient of dynamical friction, in which bmax 
enters as the argument of the logarithm function. (Note that there are other uncertainties 
as well: in the above formula for dynamical friction, the velocity dependence of the Coulomb 
logarithm was ignored by replacing it with a typical relative velocity, Vq, which is treated 
as a constant; when the velocity of the black hole is small, it is not clear that this is a 
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good approximation, and it is possible that the magnitude of the coefficient of dynamical 
friction would be modestly reduced relative to the above expressions [see Merritt 2001].) In 
this paper, we adopt a density-weighted formula for the Coulomb logarithm given by Maoz 
(1993), which provides an implicit expression for bmax] in this case, 

^-, (9) 

brmn p(0) T ' 

should replace the Coulomb logarithm in the above equations. For the case of the Plummer 
potential, we obtain 

, 1 / M \ 4 fl + 3(m/M)21 ^ ^ 

nA-^sinh-M — -- / / / ,,;J,, . 10 

\2mJ 3[l + 4(m/M)2]3/2 ^ ^ 

This gives a value for bmax that is somewhat smaller than the core radius a (a value 
often chosen for hmax)-, and indeed we find this choice provides a slightly better fit to our 
simulations (see §5) than the alternative choice of hmax ~ o- 

The third force acting on the black hole is the stochastic force denoted as F(t), which 
arises from random discrete encounters between the black hole and the stars. This force 
cannot be written down analytically in closed form, and is only defined statistically as 
described below. 

We can therefore characterize the dynamics of the black hole by 

mi^(t)+/3r(t) + A;r(t) = F(t), (11) 

which is the equation of motion of a harmonically bound Brownian particle. The spatial 
components of this linear vector equation are separable into equivalent terms, and we will 
without loss of generality concern ourselves only with its a;-component: 

mx{t) + j3x{t) + kx{t) = F^{t). (12) 

(We will use x and v^ interchangeably in the following.) This is a stochastic differential 
equation since, as noted above, the form of -Fx(^) is not known. However, since this 
stochastic force is random and rapidly varying, we expect: (1) that this force is independent 
of x; (2) that it is zero on average; and (3) that to an excellent approximation this force is 
uncorrelated with itself at different times. We may formalize these statements as 

(F,(t)) = 0, (F,(ti) F,(t2)) = C6{h - t2), (13) 

where (5 is a Dirac delta function and the angular brackets denote an average over an 
ensemble of "similarly prepared" systems of stars in each of which the black hole has 



-7- 



the same initial position and velocity. We take the factor C to be independent of r; its 
magnitude will be determined in the next section. While this definition will not allow us to 
solve equation ( P^ explicitly, we will obtain closed expressions for the time autocorrelation 
function of the black hole position and velocity in the next section. That the components of 
the random force can be separated and characterized as in the latter part of equation (|T3|) 
is at this stage an assumption; its justification must ultimately come from the agreement 
between the results of the model and the numerical simulations, as detailed in §5. 

The autocorrelation function of the stochastic force on an individual star has been 
calculated before (Chandrasekhar 1944a,b; Cohen 1975), in the approximation that the 
test star and its surrounding stars move along straight lines on deterministic orbits; in this 
approximation, the autocorrelation function falls off as slowly as the inverse of the time 
lag for a uniformly dense infinite system. However, for the case we study in this paper, 
the fall-off will be much faster as a consequence of the rapid decrease in the density of the 
system outside the core radius (see Cohen 1975), and because fluctuations will tend to 
throw the black hole and the field stars off their deterministic paths and by doing so reduce 
the correlation (see Maoz 1993). Another difference arises from the fact that we consider 
here a test object which is much more massive than the surrounding stars; since the black 
hole moves very slowly relative to the stars, in the time that the motion of the black hole 
changes appreciably, the correlations in the force due to the stars would have worn off. Our 
choice of the delta function to represent the force autocorrelation function is somewhat 
of an idealization, but is justified a posteriori by the good agreement between the model 
outlined above and the results of simulations described in §5. 

Before going on to solve the equations of motion, it is useful to list the approximations 
that have gone into setting up our model. We have assumed that a black hole of mass m is 
located near the center of a stellar system of total mass M and characteristic length scale a, 
and that the mass of individual stars m^ <^ m. Hence, the black hole's velocity is expected 
to be very small compared with the velocities of the stars, and its position is expected to 
be confined in a small region, r <^ a. We assume that the total force on the black hole is 
made up of two independent, separable parts. One (i.e., F(t)), which is due to very rapid 
fluctuations in the immediate surroundings of the black hole, is assumed to average to zero 
and to be uncorrelated with itself. The other part, which consists of dynamical friction and 
the force due to the aggregate stellar system, varies smoothly with the black hole's position 
and velocity on a time scale very much longer than that of the fluctuations. 

We assume that dynamical friction is given by the Chandrasekhar formula, which 
entails a number of additional approximations (see Tremaine & Weinberg 1984; Weinberg 
1986; Nelson & Tremaine 1999). The Chandrasekhar formula was originally derived for 



an infinite and homogeneous stellar system, but it is often employed for non-homogeneous 
systems by replacing the homogeneous density by the local density. The maximum effective 
impact parameter (for relaxation encounters between the stars and the black hole) that 
enters the Coulomb logarithm is not well-defined; we assume it to be given implicitly by the 
density-weighted expression (P) above. The gravitational encounters between the stars and 
the black hole are treated as a succession of binary encounters of short duration, i.e. as a 
Markov process. The Chandrasekhar formula approximates the orbits on which stars move 
past the black hole as Keplerian hyperbolae, even though the actual stellar orbits are more 
complex. This formula neglects the self-gravity of stars in the wake induced by the black 
hole. Despite these approximations, Chandrasekhar 's formula has been found to provide an 
accurate description of dynamical friction in a variety of astrophysical situations (see BT 
and references therein). In the present context, we will gauge the reliability of its use by 
appealing to numerical simulations to test the applicability of our model. 

We conclude this section by demonstrating that the time-scale for fluctuations in Fx{t) 
is very much smaller than the time-scale on which the position and velocity of the black 
hole change. 

Near the center of a Plummer model, where the massive black hole is localized, the 
stellar density is p ~ (SM/Ana^), since r '^ a; therefore, the typical separation between 
stars is D ~ (Ana^m^/^M) . The typical stellar velocity is Vq ~ {GM/2a) ' . The average 
time period of changes in F^fJ;), caused by discrete stellar encounters, is then approximately 
Tioc ~ D/Vq. Now the characteristic time period with which the black hole's motion changes 
is Tbh ~ Stt/cjo, where c^o = jGM/a?, as is shown in the next section. We thus have 

Tbh ^ V37r2My 1^37r2iV; A^Va' 

if there are a total of A^ stars each of mass m^ (i.e., M = m^N). Therefore, Tioc <C Tbh 
for large N, and we are justified in separating the total force on the black hole into slowly 
varying and rapidly fluctuating contributions. 



3. Solution of the Model Equations 

If we choose the initial position and velocity of the black hole to be 

a;(0) = xo, x{0) = vq, (14) 

then we can solve equation ([1^) formally as 

x{t) = e~'^^(xo cos Lut -\ sincjt) H / Fx{t — z) e~^^ sinuzdz. (15) 

^ uj '' muj Jo 
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Using Leibniz's rule to differentiate the second term on the right hand side of equation (|T5|) 
under the integral sign, we can also solve for the velocity: 

x(t) = e^'^^lvQCOsut sincijt) H / F,j.(t — z) e^'^^i—'ysinujz + ujcosuz) dz. 

(16) 
In the above equations, 

Wo = J — , 7 = 7^, ^ = \/^o - 7^- (17) 

V m Im * 

In the case of interest with ra <^ M, we have 7 <C wq, and so tu ~ t^o- Note that the exact 
results of equation (p!7| ) - and not the above approximations - have been used in comparing 
the predictions of the model with the numerical simulations of §5. 



Using the first of the properties in equation ([13D , we have the following: 



(x(t)) = xqc '''* ( cos cut H — sincijt) H e '''^sincut, (18) 

2 
(xit)) = t'oe"'^* ( cos cut sincut) — XQ—e~'^^ miujt. (19) 

Note that in the steady state (i.e., as t — > cxd), the average values of the position and 
velocity components are zero. In the above equations and the subsequent equations, angular 
brackets have the same meaning as in equation (0). 



Using the second of the properties in equation (pTSD, we can employ the delta function 
to perform the resulting double integral and solve for the time autocorrelation functions of 
the black hole position and velocity, with a time lag T: 

{x{t)x{t + T)) = {x{t)){x{t + T)) 

+ [uosinujT{l-e-'^^\l - 2 sin^ o^t + - sin 2u;t) } 



4:711?Uj'^UJn 







UJ 



2 2 

+ — coscuTil - e-2^*(l + 2^ sin^ ujt+- sin2cut)}], (20) 

7 UJ^ OJ 



{x{t)x{t + T)) = (x(t))(x(t + T)) 



+ -\- sin ujT{-1 + e"2^*(l - 2 sin^ ut - - sm2ujt)} 

2 
+ cosu;T{l - e"2^*(l + 2^ sin^ ut-- sin 2uot)}]. (21) 
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in the same way we can calculate another quantity that will be useful later: 

(x(t)x(t)) = {x{t)){x{t)) + 7r-^e-2^* ^^^2 ^^_ ^22) 

Using equations ( pOD and (pi]) , we can obtain the corresponding steady state autocorrelation 
functions for position and velocity, which are functions of the time lag T alone, by letting t 
go to infinity and thus letting transients die out. 



C^^{T) = hm (x(t)x(t + T)) = ^— ^^[cosu;T + ^ sincuT], (23) 



Ce-''^ , 



7 



C,,{T) = \im{x{t)x{t + T)) = ^[coscjT- ^sincuT]. (24) 

t-*oD A'ym'^ UJ 

Since, in most cases of interest, 7 <^ tu, these functions are essentially pure damped cosine 
terms with zero phase. Note that the stationary state autocorrelation functions are related 

by 

C....{T) = -^^^^. (25) 



It remains now to determine the constant C. If we multiply (12) by x, rearrange and 
take the ensemble average, we obtain 

{j^i^mx' + lkx')) = {F^{t)x) - iPx'). 

In the steady state, the average rate of change of total energy of the black hole is zero, and 
we get 

{F^{t)x) = {I3x^)- (26) 

in other words, the "heating" by the medium due to fluctuations must in the steady state 
equal the "cooling" due to viscous dissipation by the force of dynamical friction, which is 
a form of the general relationship between the processes of fluctuation and dissipation (see 
Bekenstein & Maoz 1992; Maoz 1993; Nelson & Tremaine 1999). 

The value of (i^) in the steady state is easily evaluated by setting T = in equation 



(H); thus. 



BT calculate the total heating per unit time to be 167r^lnA G^mm^, /^°° u/(r,M)(iti 
[adapted from equation (8-66) in BT]. Isotropy imphes that the heating due to the 
x-component alone will be a third of this quantity, namely 

16 f°° 

{Fx{t)x) = — vr^lnAG^mm^ / uf{r,u)du. 

3 Jv 
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Since the black hole velocity v is small, we can replace the lower limit in the integral above 
by zero. Then, for the Plummer model we obtain 



/ uf(r,u)duc::i Ci m \&(r) du = -Ci'^irf' c::i -Ci( for small r, 

Jv Jo \ 2/ 9 9 \ a J 



96 Ma'^ 

^1 



rVSnHGMf 

By plugging this back into the expression for (Fj.(t)i;) and using equations (P7D, ( PB| ) and 
(0), we obtain finally 

C = 'jmm^, (28) 

y (jj 

and 

,.2^ C 2GMm^ 

\X ) = = . 

47771^ 9am 

The first equality in the above equation was obtained from equation (^) by setting T = 0. 
Note that this is slightly higher than the value that would have been obtained had the 
black hole's kinetic energy been in strict equipartition with that of the stars in the core of 
the Plummer potential; had that been the case, the numerical coefficient above would have 
been 1/6 instead of 2/9 (see equation (|]), where the mean squared 3-dimensional velocity 
of the stars in the core has been calculated). 

As an aside, we point out that a similar calculation for a Maxwellian distribution of 
stars — / oc e~^ /^°" , where a is the root mean squared value of a single component of 
velocity — would have yielded 

C = Ajmm^a'^, 

and 

_ c _ ™^^ 

This is the familiar condition for equipartition between the kinetic energies of the black 
hole and a star. 

Returning to the Plummer model, we have, by making use of equations (p8[), ( |23D and 
(p^), the following expressions for the mean squared position and velocity components of 
the black hole in the steady state: 

9 m 

, „, 2GMm^ , , 

9am 
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(compare equation ( p9D to equation (101) in Bahcall & Wolf 1976; the latter is rederived in 
Lin & Tremaine 1980 just after their equation (16); these equations agree with our result 
above to within a factor of order unity). If we have A^ stars in the cluster of total mass M, 
then m^ = M/N, and we can rewrite the above equations as 

(x') - -^ (31) 

, ., 2GM^ , , 



4. The Probability Distributions of the Black Hole Position and Velocity 

Following the treatment of Chandrasekhar (1943a) and Wang and Uhlenbeck (1945), 
we can derive a partial differential equation, called the Fokker-Planck equation, for the joint 
probability distribution of the position and velocity components of the black hole. 

Let iy(x,f^,t) represent the probability distribution of the x- components of the black 
hole's position and velocity at time t; i.e., W{x, v^-, t) Ax Avx is the probability that at time 
t, the black hole lies between x and x + Ax and has a velocity between Vx and Vx + Avx- 
Let iIj^x + Ax,Vx + Avx,t + At\x,Vx,t) represent the (conditional) transition probability 
that at time t + At, the black hole is at x + Ax and Vx + Avx, given that at time t, it was 
at X and v; At is taken to be an interval that is long compared with the time-scale over 
which the stochastic force Fx{t) varies but is short compared with the time-scale on which 
the black hole's position and velocity change. 

The evolution of the probability W is expected to be governed by the following 
equation: 

W{x,Vx,t + At) = / W{x — Ax,Vx — Avx,t)^jj{x,Vx^t + At\x — Ax,Vx — Avx,t)d{Ax)d{Avx) 

Note that in writing this equation, we are assuming that the black hole's motion is a 
Markov process which depends only on its position and velocity an "instant" before, and is 
independent of its previous history. Rewriting the expression for ip in the above equation 

as ip{x — Ax + Ax, Vx — Avx + Avx, t + At\x — Ax, Vx — Avx, t), and expanding both sides 
of the equation in Taylor series, we have 

dW ^^ ^ _ f^{-l)\^ ^^ d' 

i=o ■'• 'dyi...dyi' 

where each yi is either x or f^;, and where for brevity we have defined 



W{x,Vx,t) + -^At + --- = JY.^-^Ayi...Ay,- —yj^Wix,Vx,t)]d{Ax)d{Avx 



ip^ = ip[x + Ax, Vx + Avx, t + At|x, Vx, t). 
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Keeping only derivatives up to the second order on the right hand side, we have 



dW 



ow r CI 

W{x,v^,t) + ^At + --- = J d{Ax)d{Av^)[l - Ax— - Av, 



d 



X (-\ 



The first term on the right hand side is simply W{x,Vx,t), which cancels with the same 
term on the left hand side. Dividing both sides by At and taking the limit At -^ 0, we 
obtain 

dW d d 1 d'^ 1 d'^ d'^ 

where the D coefficients are the diffusion coefficients of this general Fokker-Planck equation 
in two variables, and are defined as 



D^ = lim -— / d(Ax)d(Avr)(Ax)'d'A = lim 

At-^O At J At^O 

D,, = \im^ — J d{Ax)d{Av,){AxfiJA = jjm 
Dxv, = jim^ -r- / d{Ax)d{Avx){AxAvx)i/jA = jim^ 



(Ax) 

b At 

(AxAt;^) 
At ' 



etc. 



The diffusion coefficients can be calculated very easily by using the equation of motion 
(|l^) and the definition of the autocorrelation of the random force Fx{t) in equations ([T3|) . 
We have 

Ax = Vr^At, 

and by integrating the equation of motion for a short time At which is long enough that 
many random encounters have taken place but not so long that the black hole's x and Vx 
have changed appreciably, 

Avx = - —Vx + — X At + / ^^ dt. 
\m ml Jt m 



Based on the above and equations ([T3|) , we find that the only diffusion coefficients that do 
not vanish as At — » are: 



D. 



^X 1 ■'-^Vx 



'b k ^ 

—Vx H X 

m m J 



and D^^^^ 



C 



m 



2' 



Thus, the Fokker-Planck equation reduces to 



dW 



dx ^ dvx 



b k \ 

—Vx H X W 

m ml 



+ 



C d^W 
2w? dv^ 



(34) 
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The stationary distribution Wst{x,Vx) = W{x,Vx,t -^ oo) is found by setting the time 
derivative on the left hand side of equation ( |3lD to zero. 



The solution of equation ( plf ) is complicated, but we write it down in terms of the 
quantities derived in previous sections (see Chandrasekhar 1943a): 



W{x,x,t) 



where 



exp[—{Aw{x—{x)) + Bw (x — {x)) — 2Cw (x — (x)) {x — {x))}], 



2'K^W 



(35) 
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Note that equation (^^ describes a general Gaussian distribution in the two variables x 

and X = Vx- 

Of particular interest is the stationary distribution Wst{x.,Vx)-, which is obtained from 
equation (|35|) by taking the limit t ^ oo: 



yVst{x,Vx) = 7^— exp 
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(36) 



this is the product of two independent Gaussian distributions in the variables x and Vx- 
This is a consequence of the linear nature of equations ([TT|) and ([1^). It is easy to obtain 
the marginal stationary distributions of these variables by integrating out one or the other: 
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(37) 
(38) 



where (x^) and (f^) are given by equations (|^) and (|30D , respectively. 
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5. Tests of the Model using Numerical Simulations 

We have performed a number of computer simulations to test the validity of the model 
presented in §2-4. The code we use solves the combined dynamics of the black hole and the 
stars using different equations of motion: 

"^r = E7^r— T^TXT^w^' (39) 



k 



(|r.,fc-r|2 + £2)3/2' 



m^,kr^,k = J. 12^ 2W2 ~ "^^,fcV<l> r^,fc ; 40 

(|r-r^,fcr + e^r^ 

where m and r are the mass and position of the black hole, respectively, and m^ ^ and r^ ^ 

are the mass and position of the A;-th star, respectively; the Coulomb force is softened by 

the parameter e to prevent numerical divergences when a star passes very close to a black 

hole; and $(r) is the analytical expression for the stellar potential in equation (Q). Thus, 

the black hole interacts with the stars through a softened Coulomb force, and the stars 

interact with each other through an analytical gravitational field. Combining an analytical 

potential with the traditional "direct summation" N-body technique ensures that accuracy 

is not sacrificed in calculating the motion of the black hole (the object of greatest interest 

for us). The particles themselves are moved (with varying step-sizes which are calculated at 

every time step) using suitably modified versions of the fourth-order integrators of Aarseth 

(1994). 

The improved efficiency in the calculation is thus obtained at the price of having to 
keep the potential due to the stars (although not necessarily their density profile) fixed. 
However, this approximation does not appear to have a significant effect on our results. 
We have performed a number of simulations using other methods to test the results. 
These include the direct summation N-body code known as NBODYl (Aarseth 1994) for 
a relatively small number of particles, and the program known as SCFBDY, described 
in detail in Quinlan & Hernquist (1997). The latter program expresses the potential as 
an expansion in an appropriate set of basis functions instead of having a fixed potential 
in equation ( ^OD above; the coefficients of this expansion are self-consistently updated at 
chosen time steps. Although the precise motion of the black hole is not identical for different 
simulation methods - since the force on the black hole in each case is calculated differently 
- we have found that they all give similar results as far as the statistical properties of the 
black hole's dynamics are concerned. In particular, the mean squared values of the black 
hole's position and velocity in the stationary state of the system are approximately equal 
irrespective of the method used, and are similar to the values derived from the model 
presented in this paper. We believe that this is because the statistical properties of the 
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black hole's motion are determined primarily by the properties of the restoring force and 
dynamical friction which are provided by the unbound stars, outside the region of the black 
hole's gravitational influence. These regions are relatively unaffected by the central black 
hole if its mass is much smaller than the total mass of the stellar system. That being so, 
we have used the method of the fixed potential for the simulations described below in order 
to be able to integrate efficiently large numbers of stars for long spans of time. 

Our standard Plummer model has parameters G = M = 1 and a = 37r/16; in these 
units, the gravitational energy of the initial stellar system alone is —1/4, and the circular 
period at r = a is tdrc ~ 4.78. We take the mass of the black hole to be m = 0.01. The 
softening length was chosen to be e = 5 x 10^^. For these parameters, ujq = 2.212,7 = 0-30 
and 00 = 2.192. 

In Figure 1, we show the results of a simulation in which the black hole was started 
off with zero velocity from the origin in a system of A^ = 100, 000 stars. The first and 
second panels show the evolution of the black hole's x-component of position and velocity, 
respectively. The third panel shows the autocorrelation function of the x-component of the 
black hole's position as calculated from the simulation (the calculation was stopped at time 
t = 600), and as computed from our model; the two curves are in good agreement, at least 
for short time lags, and the discrepancies could be due to the uncertainty in the maximum 
effective impact parameter in the dynamical friction formula. Note the persistence of the 
actual autocorrelation function of the black hole, which will be discussed further below. 
The autocorrelation function of Vx is not shown; according to equation (^) it can be simply 
derived from the autocorrelation function of x by taking a double time derivative. 

In Figure 2 we test equations ( PT| ) and ( P^ which predict that the root mean squared 
position and velocity components of the black hole should decline with the total number of 
stars N as N'^/^. We show the results of 4 simulations with N = 12,500, 25,000, 50,000 
and 100,000; in each case, the simulation was stopped at time t = 600; the agreement with 
the predictions of the model is evidently good. 

In Figure 3, we test equations (0) and (pH]), which predict that the black hole's 
position and velocity components in the steady state should be Gaussian distributed; 
the empirically binned distributions were computed for the case with N = 100, 000. The 
agreement with the model predictions is again very good. 

In the above simulations, the black hole's orbit remains close to the center and appears 
to be essentially stochastic, in that it does not seem to be confined to a special sheet or 
line in phase space. At any point in time, many stars are bound to the black hole in the 
sense that they have negative energy with respect to it. Most of these stars are within the 
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gravitational sphere of influence of the black hole and their total mass is comparable to that 
of the black hole. 

The autocorrelation functions of the black hole's position and velocity do not appear 
to damp entirely with ever increasing time lag T, in contrast to equations (|23| ) and (0). 
Although Figure 1 shows only a small part of the autocorrelation function, it turns out 
that the oscillations persist for indefinitely long T at roughly the residual (and, apparently, 
somewhat varying) amplitude shown at the right of the third panel in the figure. The 
frequency of the oscillations is close to the fundamental frequency calculated in the paper 
and evident in the third panel of Figure 1. We attribute these oscillations to the presence 
of very weakly damped coherent modes in the stellar system, of the kind reported by Miller 
(1992) and Miller & Smith (1992), and calculated by Mathur (1990) and Weinberg (1994). 

While it is not the purpose of this paper to study such modes, we have, in an attempt 
to identify the source of the above oscillations, performed the following experiment. We set 
up the system of stars as above but without the black hole at its center, and kept track of 
one component of the total force at the origin of the system. The discrete Fourier transform 
of the sequence of forces at successive time steps revealed a strong peak very close to the 
fundamental angular frequency of oscillations at the bottom of the gravitational potential 



well: coq = Jk/m = JGM/a^. This could account for the undamped low-amplitude 
oscillations, at roughly the above frequency, seen in the autocorrelation functions of the 
black hole's position and velocity. Consider a simplified situation in which the black hole 
is subject to an additional force Di cos{ujot) + D2sm{uot), which is due to the conjectured 
undamped mode of this frequency mentioned above; here, Di and D2 are taken to be 
independent of time, for simplicity. If we add this force to the right hand side of equation 
(|l^), assume that it is independent of the random force Frc{t), and perform an analysis 
similar to that in §3, it is easy to see that we would obtain a new contribution to each of the 
autocorrelation functions in equations ( pB]) and ( p^ in the form of an additive term which 
is proportional to cos{uJoT), i.e., a term that does not damp with increasing T (note that 



for most systems, uq is approximately equal to u; = Juq — 7^, the frequency with which 
the black hole's position and velocity autocorrelation functions oscillate). The addition of 
such a term would not affect the good agreement for small time lags between equations 
(^) and (p4[) and the results of numerical simulation (Figure 1), since the amplitude of 
these residual oscillations is very much smaller than the amplitude of the autocorrelation 
functions for small T. 
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6. Lower Limit on the Mass of Sgr A* 

We may use our model to derive a lower limit on the mass of the black hole in the 
Galactic center, Sgr A*. The observed upper limit of 20 km s~^ on the intrinsic proper 
motion of Sgr A* (Reid et al. 1999), when combined with equation (|30|) , provides such a 
limit. 

Measurement of proper motions of stars close to Sgr A* indicate that a total mass of 
~ 2.6 X 10^ Mq resides within a distance of ~ 0.01 parsec of Sgr A* (see, e.g., Eckart & 
Genzel 1997, Ghez et al. 1998, Ghez et al. 2000). Not all of this mass need be attributed 
to the black hole; some of it could be due to a cluster of stars surrounding it. If we assume 
that this cluster is distributed according to a Plummer profile, then we can apply the results 
derived in previous sections to the entire system comprised of the stellar cluster and the 
black hole. 

Let us set the total mass inside a distance R = 0.01 pc from Sgr A* to be 
Mn = 2.6 X IO^Mq. Using m as the mass of the black hole, we have the condition 

where the second term is the mass of the stellar cluster (of total mass M) within R. 
Combining this with equation (|30|) in the form 



we obtain 



9am 



m > —^ j-^, (42) 

^ ^ 2 G(i?2+a2)3/2m^ 

where < f ^ >rnax is the maximum mean squared speed of one component of the black hole's 
velocity. We can now set {v'^)Uax ~ 20 km s~^ to get an approximate lower limit on m, 
assuming that m^ ~ IMq. This relation is plotted in Figure 4(a) for various values of a, the 
scale length parameter of the Plummer cluster. The mass of the black hole must be given 
by points lying above the curved solid line. Evidently, this relation implies a lower limit on 
m. 

We can derive stricter limits on m by noting, as did Maoz (1998), that the allowed 
values of a are restricted by the condition that the upper limit of the lifetime of a cluster 
of stars is its evaporation time, when stars would have escaped from the cluster because 
of scattering. The evaporation timescale of the cluster should be long enough to make it 
probable that the cluster be observed at the present epoch. A reasonable assumption is 
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that this timescale, t^vap, is bounded between the values 1 Gyr and 10 Gyr, the latter being 
the approximate age of the Galaxy. The evaporation timescale is tevap ~ 300tr/i, where trh 
is the median relaxation time, given by (see, e.g., Maoz 1998, BT) 



0.14A^ 



\n{0 AN) V CM' 

where A^ = M/m^, is the number of stars in the cluster (with M given by equation ^l]), and 
Th is the system's median radius; Vh = 1.3a for the Plummer profile. 

Setting tevap = 1 Gyr, we obtain the dashed line in Figure 4(a); tevap > 1 Gyr denotes 
the region to the right of this line. Hence, the allowed values of m and a lie to the right of 
this dashed line and above the solid line. The minimum value of m under these assumptions 
is given by point A at which m ^ 3.8 x IO^Mq, a ~ 0.021 pc, and the total mass of the star 
cluster is M ^ 3.3 x IO'^Mq. 

If we perform a similar calculation for tevap > 10 Gyr, we obtain the dotted line in Figure 
4(a); the minimum value of m is then given by point B at which m ~ 1.6 x W^Mq, a ~ 0.049 
pc and M ^ 3.3 X IO^Mq. 

Reid et al. (1999) expect further observations to reduce the limit on the peculiar motion 
of Sgr A* from 20 km s^^ to 2 km s~^ over the next few years. Figure 4(b) repeats the 
above calculations for this limit. The condition tevap > 1 Gyr then gives a new lower limit 
for m (point A), at which m ^ 3.4 x IO^Mq, a ^ 0.022 pc, and M ^ 3.0 x WMq. Point B, 
the minimum of m if tevap > 10 Gyr, is now characterized hj m ^ 1.1 x IO^Mq, a ~ 0.054 
pc, and M ^ 2.5 x IO^Mq. In the latter case, the attainable lower limit on m will be 
interestingly close to its upper limit of Mfj = 2.6 x 1O^M0. 



7. Summary 

In this paper, we have developed a stochastic model to describe the dynamics of a 
black hole near the center of a dense stellar system. The total force on the black hole is 
decomposed into a slowly varying part originating from the response of the whole stellar 
system, and a random, rapidly fluctuating part originating from discrete encounters with 
individual stars. We have shown that the time scale over which the latter force fluctuates 
is very short compared with the time scale over which the former changes; hence the 
justification for the separation of the total force into these two independent components. 
The slowly varying force itself is approximated as the sum of two contributions: the force 
on the black hole due to the potential of the whole stellar system and a force of dynamical 
friction which causes the black hole to decelerate as it moves through the stellar system. 
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The stochastic force at the position of the black hole is assumed to have a zero average 
and to be essentially uncorrelated with itself over time scales that are short compared with 
the characteristic period over which the velocity of the black hole changes considerably, but 
long enough for many independent fluctuations of the stochastic force to have occurred. 

If the stellar system is approximated by a Plummer model, then the problem essentially 
reduces to describing the Brownian motion of a particle in a harmonic potential. We have 
shown that after long times, the black hole's velocity has a zero average and its average 
location coincides with the center of the stellar potential. However, the root mean squared 
position and velocity tend towards non-zero values which are independent of time and 
of the black hole's initial position and velocity. The steady-state time autocorrelation 
functions of the position and velocity were shown to be approximately given by damped 
cosine functions. For a Maxwellian distribution of stars, strict equipartition of kinetic 
energy between the black hole and the stars is achieved. 

The model was completed by solving for the one unknown parameter — the mean 
squared magnitude of the stochastic force on the black hole — by making use of the 
close relationship between processes of fluctuation and dissipation, according to which the 
"heating" by the fluctuating force equals the "cooling" caused by the dissipative force of 
dynamical friction (the third force on the black hole — the force due to the stellar potential 
— is conservative). 

A Fokker-Planck equation for the diffusion of the probability distribution of the black 
hole's position and velocity was developed and solved. The solution implies that in the 
steady state these variables are distributed independently as Gaussians. 

The predictions of the model were tested by comparing with the results of various 
N-body simulations; the agreement is good, thus justifying the elements of the model. 
In the simulations, we find possible signs of the existence of weakly damped coherent 
modes associated with the stellar system (Weinberg 1994). The total force at the origin 
of the system has a strong Fourier component at the frequency of fundamental oscillations 
characteristic of the approximately harmonic form of the stellar potential near the center. 
If we consider the black hole to be a test particle affected by such a force, then the 
autocorrelation functions of the black hole's position and velocity would be found not 
to damp with the time lag, but to persist; the results of simulation do indeed show the 
presence of such oscillations at very low amplitude for arbitrarily long time lags. 

Finally, we applied the results of the model to Sgr A*. Observational limits on 
the peculiar motion of Sgr A* were used to obtain a lower limit on its mass, under the 
assumption that it is localized near the center of a system of equal mass stars distributed 
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according to the Plummer model. More stringent limits were then deduced by requiring 
that the evaporation timescale of the cluster of stars be larger than 1 Gyr. 

The Plummer model is a reasonable choice for the black hole problem considered here 
since it results in a separable system of linear stochastic differential equations; it is not clear 
that arbitrary potential-density pairs would give equations that are similarly tractable. In 
particular, what would happen in the case of black holes at the centers of galaxies if those 
are described by singular power-law density profiles? A simple-minded generalization of 
our method would not necessarily work (for example, because of possible divergences in the 
distribution function). However, we believe that our model captures the qualitative features 
of more complicated situations for two related reasons. First, we note that the Plummer 
model is not an equilibrium solution for a cluster of stars with a black hole present (see, e.g., 
Huntley & Saslaw 1975 and Saslaw 1985); even in non-singular models such as the Plummer 
model, the black hole ultimately induces a density cusp. We have ignored this complication 
in our model and found that the model nevertheless provides a good description of detailed 
numerical simulations. Second, the black hole tends to carry its cusp of bound stars with it 
as it moves around; thus, it is as if a black hole of a somewhat larger effective mass were 
moving in a background consisting of unbound stars. With the cusp effectively removed, 
the density profile of this background would be flat near the center. Since the restoring 
force and dynamical friction are provided mainly by the unbound stars, we believe that the 
essential components of our model are still valid. (For similar reasons - see §5 for details 
- the fact that our simulations use a fixed stellar potential is not expected to alter our 
conclusions.) We have carried out numerical simulations for the case of a particular density 
profile with a singularity, namely the Hernquist (1990) model. We find that the qualitative 
results are similar to those described here. In particular, the early-time autocorrelation 
functions of the black hole's force and velocity continue to be described well by damped 
cosine functions of fixed frequency. The detailed characterization of the black hole behavior 
in terms of the model parameters is different, and requires a more careful calculation. 
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Fig. 1. — Results from a simulation of a black hole of mass m = 0.01 in a Plummer model 
of total mass M = 1 and N = 100, 000 stars. The evolution of the x-components of the 
black hole's position and velocity are shown from t = to t = 100. The third panel shows 
the autocorrelation function of x for time lags from T = to T = 20. The solid line is the 
autocorrelation function obtained from the simulation, and the broken line is that predicted 
by our model. 
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Fig. 2. — Results from simulations of a black hole of mass m = 0.01 in a Plummer model 
of total mass M = 1 and iV = 12,500, 25,000, 50,000 and 100,000 stars. The panels show 
the root mean squared values of the x-components of the black hole's position and velocity. 
The crosses show the results obtained from the simulations, and the solid lines are the values 
predicted by our model. In each case, the simulation was stopped at time t = 600. 
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Fig. 3. — Results from a simulation of a black hole of mass m = 0.01 in a Plummer model 
of total mass M = 1 and A^ = 100, 000 stars. The distributions of the x-components of 
the black hole's position and velocity are shown as histograms. The solid lines show the 
corresponding bin values according to the analytic model of this paper. 
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Fig. 4. — The region above the sohd hne in the upper panel denotes the allowed values of 
m, the mass of the putative black hole in the Galactic center Sgr A*, and a, the length scale 
characterizing the cluster of stars (each of mass 1 Mq) surrounding Sgr A*, given that the 
rms value of one component of the black hole's velocity is < 20 km s~^. The region to the 
right of the dashed line denotes the allowed values of m and a if the evaporation timescale 
of the cluster, tevap, exceeds 1 Gyr, and the region to the right of the dotted line denotes 
the allowed values if tevap > 10 Gyr. The lower panel is the same as the upper panel, except 
that the rms value of one component of the black hole's velocity is now assumed to be < 2 
km s^^. 



